library(Biostrings)
library(XML)
library(RMySQL)
library(GenomicRanges)
library(ggbio)

baseDir <- "/home/haseong/alpha/dev/FT-NGS"
targetDir <- "/home/haseong/alpha/dev/FT-NGS"
dbName <- "/home/haseong/alpha/apps/blast/db/nt"
annodb.dir <- "/home/haseong/alpha/apps/blast/microannodb/"
setwd(baseDir)

file.names<-dir(pattern="^KHS.*fna$")
dir.names<-dir(pattern="^out_KHS.*fna$")

##for(fn in file.names){
fn<-file.names[2]
  contig.file<-paste(baseDir, "/", fn, sep="")
  contig.seqs<-readDNAStringSet(contig.file, "fasta")
  contig.names<-unlist(lapply(strsplit(names(contig.seqs), split="\\s"), function(x){return(x[1])}))
  names(contig.seqs)<-contig.names
  out.dir<-paste(targetDir, "/out_", fn, sep="")
  summary.file.names<-dir(out.dir, pattern="^con.*RData")
  
  contig.cnt<-0
  for(sf in summary.file.names){
    sf<-summary.file.names[1]
    load(paste(out.dir, "/", sf, sep=""))
    sel.contig.name<-unlist(strsplit(sf, split="\\."))[1]
    gi.list<-unlist(lapply(strsplit(as.character(bl[,3]), split="\\|"), function(x){return(x[2])}))
    gis<-names(table(gi.list))
    
    anno.matrix<-matrix("", 5000, 7)
    colnames(anno.matrix)<-c("gi", "gene", "locus_tag", "product", "protein_id", "start", "end")
    contig.cnt<-contig.cnt+1
    gi.cnt<-0
    cnt<-1
    for(sel.gis in gis){
      gi.cnt<-gi.cnt+1
      ##sel.gis<-gis[2]
      ##if(is.na(match(dir(annodb.dir), sel.gis))){
      cat("Download annotation info... GI:", sel.gis, "(", gi.cnt, "/", length(gis), "summary:", sf, "(", contig.cnt, "/", length(summary.file.names), ")", "file:", fn, "(", length(file.names), ")", "\n");flush.console()
      ##run.string<-paste("curl -s \"http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=", sel.gis, "&rettype=gb\" | awk \'/ CDS/,/db_xref=/\' | sed \'s/^\\s*//\' > ", annodb.dir, sel.gis, sep="")
      run.string<-paste("curl -s \"http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=", sel.gis, "&rettype=gb\" | awk \'/( CDS|LOCUS)/,/db_xref=/\' | sed \'s/^\\s*//\'", sep="")
      tmpr<-system(run.string, intern=T)  
      ##cat("Update mysql... ");flush.console()
      ##tmp.read<-readLines(paste(annodb.dir, sel.gis, sep=""), n = -1)
      cur<-0
      for(i in 1:length(tmpr)){
        sel.line<-gsub("\"", "", tmpr[i])
        if(cur==0){
          if(length(grep("CDS", sel.line))>0){
            m<-gregexpr("[0-9]+", sel.line)
            pos<-as.numeric(unlist(regmatches(sel.line, m)))
            anno.matrix[cnt,"start"]<-min(pos)
            anno.matrix[cnt,"end"]<-max(pos)
            cur<-1
          }else if(length(grep("LOCUS", sel.line))>0){
            tmpstr<-unlist(strsplit(sel.line, split=" "))
            ref.seq.len<-as.numeric(tmpstr[which(tmpstr=="bp")-1])
          }else if(length(grep("ORGANISM", sel.line))>0){
            organism.name<-gsub("^\\s*", "", gsub("ORGANISM", "", sel.line))
          }
          
        }else{
          if(length(grep("/gene=", sel.line))>0){ 
            anno.matrix[cnt, "gene"]<-unlist(strsplit(sel.line, split="gene="))[2]
          }
          if(length(grep("/locus_tag=", sel.line))>0){ 
            anno.matrix[cnt, "locus_tag"]<-unlist(strsplit(sel.line, split="locus_tag="))[2]
          }
          if(length(grep("/product=", sel.line))>0){ 
            anno.matrix[cnt, "product"]<-unlist(strsplit(sel.line, split="product="))[2]
          }
          if(length(grep("/protein_id=", sel.line))>0){ 
            anno.matrix[cnt, "protein_id"]<-unlist(strsplit(sel.line, split="protein_id="))[2]
          }
          if(length(grep("/db_xref=", sel.line))>0){ 
            ##anno.matrix[cnt, "contigs"]<-sf
            anno.matrix[cnt, "gi"]<-sel.gis
            if(cnt%%5000==0){
              anno.matrix<-rbind(anno.matrix,  anno.matrix<-matrix("", 5000, 7))
            }
            cnt<-cnt+1
            cur<-0
          }
        }
      }
      anno.matrix<-anno.matrix[anno.matrix[,1]!="",]
    }
    
    save(anno.matrix, file=paste(out.dir, "/anno_", sf, sep=""))
  }
  
##}